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Abstract 

Recent experimental studies of molten metal droplets wetting high temperature reactive sub- 
strates have established that the majority of triple-line motion occurs when inertial effects are 
dominant. In light of these studies, this paper investigates wetting and spreading on reactive 
substrates when inertial effects are dominant using a thermodynamically derived, diffuse interface 
model of a binary, three-phase material. The liquid-vapor transition is modeled using a van der 
Waals diffuse interface approach, while the solid-fluid transition is modeled using a phase field 
approach. The results from the simulations demonstrate an O (t~^ 2 ^ spreading rate during the 
inertial regime and oscillations in the triple-line position when the metal droplet transitions from 
inertial to diffusive spreading. It is found that the spreading extent is reduced by enhancing 
dissolution by manipulating the initial liquid composition. The results from the model exhibit 
good qualitative and quantitative agreement with a number of recent experimental studies of high- 
temperature droplet spreading, particularly experiments of copper droplets spreading on silicon 
substrates. Analysis of the numerical data from the model suggests that the extent and rate of 
spreading is regulated by the spreading coefficient calculated from a force balance based on a plausi- 
ble definition of the instantaneous interface energies. A number of contemporary publications have 
discussed the likely dissipation mechanism in spreading droplets. Thus, we examine the dissipation 
mechanism using the entropy-production field and determine that dissipation primarily occurs in 
the locality of the triple-line region during the inertial stage, but extends along the solid-liquid 
interface region during the diffusive stage. 
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I. INTRODUCTION 



Characterizations of metal alloys wetting and spreading on dissolving substrates typi- 
cally assume that inertial effects are not dominant or that the majority of dissipation is due 
to viscous forces [IHI]. In many respects this seems an entirely reasonable approach since 
the majority of experiments do not capture the early time behavior when inertial effects 
are dominant, but focus on the late-stage spreading when chemical-diffusion dominates and 
substrate dissolution occurs. Typically, experimental studies measure only slow spreading 
on the order of seconds or even minutes for millimeter-sized metal droplets consistent with 
diffusion-dominated spreading [TJ |SH7j. However, using improved techniques, a number of 
recent experiments [7HTU] capture the rapid early-stage spreading and demonstrate that the 
spreading duration is consistent with the inertial time scale pjQ. The variations in exper- 
imental findings can be attributed to differences in substrate temperature, composition of 
the vapor phase influencing substrate oxidation, contact mechanisms between the substrate 
and molten droplet, camera shutter speed, as well as other factors [7j. An often important 
aspect of managing these factors is arresting the formation of a substrate ridge on which the 
triple line becomes attached, which can retard spreading considerably [5]. 

The spreading droplet is often characterized in terms of a velocity versus contact angle 
relationship where the velocity is scaled using the instantaneous Capillary number, Ca* = 
U*is/j, where U* is the instantaneous spreading speed, v is the liquid viscosity and 7 is 
the liquid-vapor interface energy. Saiz et al. postulated that the dissipation mechanism 
may not be due to viscous forces as previously understood [3, [12] . Clearly, in cases where 
the dissipation mechanism is not due to viscous effects, Ca is no longer a useful quantity 
for characterizing the spreading and an alternative parameter is required. An effective 
"triple-line friction" derived from molecular kinetics theory is suggested by Saiz et al. that 
is independent of viscosity but still dependent on interface energy and the contact angle. 
A number of recent experimental studies [TT] clearly show that a large proportion of the 
spreading is characterized entirely by the inertial time scale (U = \J pB^/'y, where p is the 
liquid density and Ro is the drop radius) with U ~ t~ 1//2 , which is much faster than typical 
viscous spreading laws [13J. Furthermore, molecular dynamics studies of Ag-Ni and Ag-Cu 
systems seem to confirm the t~^ 2 dependence of the spreading rate even for relatively small 
droplets HUE*]. 
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This paper employs a diffuse interface method in order to analyze the issues surround- 
ing the inertial spreading regime and dissipation mechanism discussed above. The diffuse 
interface approach implicitly includes a wide range of phenomena and as such does not re- 
quire a posited relationship between spreading rate and contact angle [16J. Villanueva et 
al. jl] used a diffuse interface method to model reactive wetting and clearly identified two 
separate spreading regimes: an initial viscous regime and a subsequent diffusive regime [2]. 
The viscous regime demonstrated excellent agreement with standard viscous spreading laws. 
Further work by these authors [I] employed the same model to examine the effects of dis- 
solution on spreading by first recovering the non-dissolutive hydrodynamic limit as a base 
state. In the viscous regime they found the spreading to be independent of the diffusion 
coefficient, but accelerated in the diffusive regime as the diffusion coefficient is increased. 
This paper outlines a similar process using the initial liquid concentration to vary the driving 
force for dissolution, while maintaining a constant diffusion coefficient. The general consen- 
sus of the literature is that inertial spreading occurs more slowly in systems that exhibit 
dissolution than in immiscible systems that do not exhibit dissolution [Tj [T7] . However, this 
is contradicted by a number of experiments for saturated and pure liquids that show that 
the spreading can be on a similar time scale under certain experimental conditions [TJ [9] . 

The work of Villanueva et al. jlj considers droplets that do not exhibit inertial effects 
due to the small drop size, which is limited by the requirement of having a narrow interface 
(« 1 nm). In contrast to reference [I], this work sacrifices the realistic interface width in an 
attempt to model a system that exhibits inertial effects. Due to the drop size restrictions, 
the inertial time scale used in Villanueva et al. is ti ~ 6 x 10~ n s and the capillary time 
scale, t c = uRq/j 2 x 1CT 11 s. At these values, the extent of spreading during the inertial 
stage is limited and the characteristic inertial effects are suppressed by viscous forces. The 
Ohnesorge number, given by Oh = t c /ti, quantifies the relative importance of inertial and 
viscous effects. Typically, millimeter-sized metal droplets are highly inertial in nature with 
Oh pa 1 x 1CT 3 . Characteristic inertial effects, such as triple-line position oscillations and 
large droplet curvature variations, are reduced for Oh > 0.01 and eliminated for Oh > 1 [18J. 
In Villanueva et al, Oh w 0.3 and in this work Oh w 6 x 10~ 3 . 

Jacqmin makes an extensive study of the role of the diffuse interface method, specifically 
for a Cahn-Hilliard-van der Waals system (CHW), in relieving the stress singularity that 
occurs for classical sharp interface methods [TB]. Since the interface is diffuse, the CHW does 
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not require an explicit alteration to the no-slip boundary condition to allow for triple-line 
slip. Jacqmin demonstrates that the CHW has the same far field and macroscopic behavior 
as classical hydrodynamic models of slip. Thus, in diffuse interface models that include 
hydrodynamics there is no need to define a slip length. The interface width determines both 
an effective slip length and the concentration profiles within the diffuse interface associated 
with adjustments to adsorption and desorption; these factors affect the evolution of the 
system in subtle ways. There is no exact expression relating interface width and the effective 
slip length, however, A = 5/2R is suggested as a good rule of thumb in Ding and Spelt [T9"] . 
where A is the dimensionless effective slip length for a diffuse interface model. It is claimed 
that the slip length can be as large as 50 nm [20] , which is close to the chosen interface width 
in the present work, although the drop radius is only 1 ]±m. The slip length is found by Ding 
and Spelt to influence the onset of oscillations that occur when the droplet transitions from 
the inertial stage to the diffusive stage. The critical value of Re* for which oscillations occur 
is reduced with decreasing A. Hocking and Davis [21] have demonstrated that there is no 
simple relationship between contact angle and velocity when the approach to equilibrium 
becomes oscillatory, which seems to be the case in a number of experimental and numerical 
studies of millimeter-sized droplets [8] [91 [18], [19] . 

The code used for the numerical analysis in this paper is developed using the FiPy PDE 
solver [22]. Details of how to install FiPy as well as the reactive wetting code used here are 
given on the FiPy web site [23] • The numerical analysis and figures presented in this paper 
can be reproduced with the open source tools available. The underlying linear solvers and 
parallel capabilities are provided by the Trilinos tool suite [21] . 

In the following section the governing equations are presented followed by a discussion of 



the associated dimensionless parameters in section III Results from the numerical solution of 



the governing equations outlined are presented in section [TV] Section [V] analyzes the results 
in the context of previous work and ends with a discussion of the dissipation mechanism. 
Section VI presents the conclusions. Appendix |A| derives the governing equations presented 
in section [TTJ while appendix [B] presents details of the numerical methods. 
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II. GOVERNING EQUATIONS 



In this section, the final forms of the governing equations are presented along with the 
associated thermodynamic parameters and functions. The full derivation of the governing 
equations is described in appendix [Aj The system consists of a three phase (solid, liquid 
and vapor) binary alloy. The liquid-vapor system is modeled as a two component van der 
Waals fluid, while the solid-fluid system is modeled with a phase field description. The 
density field acts as the order parameter for the liquid-vapor transition. Thus, the system 
is fully characterized by the spatio-temporal evolution of the mass density of component 
1, pi, the mass density of component 2, p 2 , the phase field, 0, as well as the barycentric 
velocity field u, as determined through the momentum equation. The three dimensional 
equations are reduced to two dimensions by imposing cylindrical symmetry about r = 0. 
The initial configuration consists of a spherical droplet with a radius of 1 /im tangent to 
a solid substrate surrounded by a vapor. The incompressible approximation is not made 
in this work for numerical reasons outlined in appendix |Bj all the phases are compressible. 
The solid is modeled as a very viscous fluid as in previous phase field reactive wetting 
studies [2111]. As the total mass density, p = pi + P2, appears so frequently in the equations, 
it is more convenient to use p and p 2 as the independent density variables. For economy in 
notation, we write spatial derivatives di = d/dxi, df = d 2 /dxf and require that repeated 
indices are summed, unless otherwise indicated. Note that although the equations are solved 
with cylindrical symmetry, the equations are presented in the following Cartesian forms: 

1. Continuity 




(1) 



2. Diffusion 
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3. Phase 



— + = e^M^. (j) ~~rp~ ~QA ) ( 3 ) 



^. Momentum 



+ dj (pu iUj ) = dj (u [djm + d iUj ]) 



(4) 



where is a velocity component, T is the temperature and M = Mpip 2 /p 2 is the chemical 



mobility, which is proportional to the diffusivity, D, as outlined in Eq. (14). The values of 



M and v vary from the solid to the fluid phases with the interpolation scheme chosen to be 

M = MfM)-* (5) 

and 

v = ufu}^ (6) 

where if) = <p a with a = 4. The values used in the simulations for M s , Mf, v s and Vf are in 
Table [Tj The choice of a is discussed in subsection V C The free energy per unit volume is 
postulated to have the form 



/ = p(4>) fs + (i - p{4>)) f f + w<p 2 (l - 

where W is the phase field barrier height and p(0) = 3 (1O— 150+60 2 ) represents a smoothed 
step function common in phase field models [26]. The free energies per unit volume in the 
separate fluid and solid phases are given by, 

, eip? e l2 pip 2 e 2 p 2 2 RT 

// = — 7r H 5 1 it H [pi mpi + p 2 lnp 2 - pm(m - vp)\ (7) 

m Z m Z m A m 



and 



A x p x A 2 p 2 RT B , ref . 2 

/ s = + + pilnpi +p 2 lnp 2 -plnp + [p s - p) 8 

m m m pm v 7 
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where m is the molecular weight (assumed to be equal for each component), R is the gas 
constant, v is the exclusion volume due to the finite size of the atoms, B is the solid com- 
pressibility, pi is a reference density for the solid and the CiPi/m are the free energy con- 
tributions per unit mole due to intermolecular attraction in the van der Waals model. The 
A\ and A 2 are temperature dependent parameters related to the heat of fusion between the 
solid and fluid phases. Along with the free energy, the specification of the pressure and the 
non-classical chemical potentials are required to fully define the system, 

tf C = |£ - ^d] Pl (10) 

»Z C = |^-^j P2 (ii) 

OP2 

where €\ and e 2 are free energy gradient coefficients. The parameter values for Eqs. Q 
to (11) are presented in Table [Tj The corresponding isothermal phase diagram for the molar 



fraction of component 1 verses the molar volume is displayed in figure [T} 

Eqs. Q-Q are solved using a cell-centered, collocated finite-volume (FV) scheme. The 
solution algorithm uses a fully coupled Krylov solver with Picard non-linear updates using 
Pi, P2, <f) an d u as the independent variables. Further discussion of the numerical approach 
is given in appendix [Bj 

III. DIMENSIONLESS EQUATIONS AND TIMESCALES 

It is useful for the purposes of analysis and completeness to clearly present the various 
dimensionless numbers and time scales that arise from solving Eqs. (fi]) ([2| ^ and ^ in 
the context of spreading droplets. The dimensionless forms of Eqs. ^ and ^ are given by, 



°§ + dj (u jP2 ) = ±d s {^0, (p 2 - p, -Qdl (p 2 - Pl ))) (12) 



and 



(13) 

where the variables and operators are now dimensionless (the analysis of Eqs. ([T]) and ^ 
is not particularly revealing and is omitted). For completeness, all the time scales referred 



to in this paper are displayed in table [IT] as a prerequisite for presenting the dimensionless 



numbers in table III It should be noted that in table |H] U* — U* (t) is the instantaneous 
spreading speed and U is a fixed spreading speed posited a priori. 

The time scale idiff represents the time required for the solid-liquid interface to move a 
distance 5 due to diffusion mediated melting or freezing. The expression for t<as = 5 2 /AK 2 Df 
is determined using an error function based similarity solution (see (27]) where K is the 
solution to 

K | ( X{-X l r« \ exp(-if 2 ) 1 
+ ^ X{-X[ J 1 -erf (#) 

and the chemical diffusion coefficient in the fluid, Df, is defined by 

D f = = 9.58 x HT 10 m 2 /s (14) 

If we substitute Ro for S in the expression for idiff, a rough estimate is obtained for complete 
equilibration of the system. Since t^is ^> U, the motion of the solid-liquid interface is 
negligible for a simulation that is both computationally feasible and adequately resolves the 
inertial time scale. The motion of the solid interface due to dissolution is controlled by both 
diffusion (tdiff) and boundary kinetics (represented by t^). Here <C t i: thus dissolution will 
be limited by diffusion rather than boundary kinetics. Additionally, solid interface motion 
due to hydrodynamic effects is negligible because the solid viscosity is chosen such that 
t s 3> ti where t s represents the time scale for discernible motion of the solid. 



Table |III| presents the dimensionless numbers in terms of their constituent time scales 
where appropriate. Note that there are now two separate expressions for both the Reynolds 
number and the Capillary number based on U and U*. By making an informed choice for 
the value of U, estimates are obtained for the likely values of the dimensionless numbers 
when using U*. Here, U = Ro/ti = 5.08xl0 1 m/sis selected based on the spreading 
rate for a system that is dominated by inertial effects. The values of Oh, Re and Pe in 



table III all indicate that the interface energy and inertial forces dominate over viscous and 
diffusive forces. Since We = 1, the interface energy and inertial forces are of approximately 
equivalent magnitude. Small values of Oh are representative of many experimental systems 
of technical interest: for example, Oh f» 2 x 10~ 3 for a millimeter sized droplet of copper 
and Oh ^2x 10~ 2 for a micrometer-sized drop of lead. 



Other dimensionless numbers (included for completeness) in table III include the Mach 
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number, Ma, which requires a definition for the speed of sound in the liquid, given by [28J, 

= 8.89 x 10 2 m/s 




e< 1 1 1 

Pi 



and Q, which represents the ratio between interface and internal forces in the liquid droplet, 
but has not been identified in the literature by the authors. 

IV. RESULTS 

In this section, we explore the rate and extent of droplet spreading based on variations in 
the initial liquid concentration and the Ohnesorge number. The initial liquid concentration 
determines the driving force for dissolution, while manipulating the Ohnesorge number in- 
fluences the impact of inertial effects on spreading. The results presented here will provide 
the basis for comparison with other authors' work in section |V} 

The extent of dissolution is established by decreasing the initial value of the liquid con- 
centration, X[, requiring the solid to dissolve in order to restore X[ to its equilibrium value, 
X[' equ . Explicitly, wc set 

x[ (t = o) = (i - o ^ equ 
Pl (t = o) = P r 

where £ defines a measure of the magnitude of the driving force for dissolution (£ < 
induces freezing). When £ = 0, the system has no potential for dissolution, similar to pure 
hydrodynamic spreading where surface tension forcing dominates and interface motion is 
due only to convection as phase change is negligible. In this limit, comparisons can be 
made with simpler spreading models and power laws. In addition to the hydrodynamic case 
(£ = 0), simulations were conducted with values of £ = 0.5 and £ = 0.9. 

Figure [2] demonstrates the highly inertial nature of the spreading dynamics. Upon initi- 
ation of the simulation, pressure waves appear at the interface regions and travel through 
the interior of the droplet, but then disperse quickly. Simultaneously, triple-line motion 
begins with a rapid change in the local contact angle, but without any discernible motion 
elsewhere on the drop interface. This initiates the most noticeable feature of the spreading: 
a capillary wave propagates from the triple line along the liquid-vapor interface, initiating 
the onset of the triple-line motion and progressing to the top of the droplet, causing a rapid 
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rise in the drop height. The wave then travels back to the triple-line location while the 
droplet completes the majority of the spreading, with both events having a duration that 
corresponds to ~ 2tj. During this interval, the triple-line motion is monotonic and with- 
out interruption. On return to the triple-line location, the wave induces a reversal in the 
triple-line motion. Subsequent waves induce further reversals in the triple-line motion and 
the drop height with a period of ~ 2ij. The amplitude of the oscillations diminishes in the 
manner of an under-damped oscillator, completing approximately 5 or 6 full cycles before 
ceasing entirely. Subsequently, very slow monotonic spreading occurs with the liquid-vapor 
interface appearing to have almost constant curvature. 

Figures [3] and |4] display the scaled radial position of the triple-line, ru/Ro, against the 
scaled time for varying values of £ and Oh. The two intervals of fast and slow monotonic 
spreading can clearly be seen as well as the intervening period of oscillatory spreading as 
discussed in the previous paragraph. Increasing £ reduces the extent of spreading slightly, 
while increasing Oh eliminates the oscillations entirely and considerably reduces the spread- 
ing rate. In each of these cases, the actual amount of substrate dissolution is negligible (the 
solid- fluid interface moves less than 5/5) due to the large disparity between the dissolutive 



and inertial time scales as discussed in section III In figure [3j at early times (t < O.ltj), 
the value of £ has no impact on the spreading, but at later times (t > O.ltj) the curves 
diverge. When t > 10ij, the curves stop diverging and seem to remain at a fixed distance 
apart. Increasing £ not only results in a slight reduction in the extent of spreading, but also 
results in a notable reduction in the amplitude of the oscillations. These factors indicate 
that there is a seemingly modest decrease in the driving force for spreading with increasing 
£. In figure [4j the Oh = 5.7 x 10 _1 curve diverges from the other curves at very early times 
and has a greatly diminished spreading rate. Eventually, the curves become coincident at 
late times when the spreading is free of observable inertial manifestations for all values of 
Oh. 

In order to compare with other models, the radial position results presented in figure [4] 
are presented using a scaled spreading velocity in figure [5j The spreading velocity is scaled 
using a Reynolds number, A Re* (A = 5/2Rq), based on the interface width, 5, rather than 
using a standard Reynolds number based on the initial drop radius, Rq [T9]. The spreading 
velocity data used in figure [5] is smoothed to remove noise on the order of a grid spacing, 
the details of which are described in appendix [Bj The sign changes in the blue curve, when 
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ti < t < 10£j, correspond to the triple-line oscillations seen in figure [4| The oscillations lie 
between intervals with monotonically decreasing spreading velocity. The Oh = 5.7 x 1CT 1 
(yellow) curve exhibits a fairly steady decrease in velocity and then a much sharper reduction 
when t ~ 10tj, which corresponds to a slope change in frigure |4j Note that the Oh values 
for simulations presented in figure [5] are manipulated by changing the value of Uf only, and 
thus, a corresponding figure with no scaling for the spreading velocity would show only slight 
differences between the vertical positions of the curves. 

Figure [6] displays the apparent contact angle, 9, against the Capillary number for Oh = 
5.7 x 1CT 3 and £ = demonstrating the convergence of 9 to the nominal flat-interface, 
equilibrium contact angle, 9 eqn . The angle 9 is calculated using techniques similar to those 
described in Villanueva et al. jl]. Although 9 exhibits a hysteresis loop, it remains relatively 
steady during the period of oscillatory spreading and only varies by ~ 0.037T rad for the 
largest oscillation. 

V. DISCUSSION 

A. Comparison with other models 

At early times, the flow is dominated by inertia and comparisons with theories of spread- 
ing on flat, non-reactive substrates are fruitful. Indeed, an analytical spreading rate for 
the inertial regime can be derived, see Biance et al. [13], and is given by i^ 1 / 2 . In fig- 
ure [5j the slope of this power law (black dashed line) shows reasonable agreement with the 
Oh = 5.7 x 10~ 3 (blue curve) during the inertial regime. The vertical position of the black 
dashed line is selected to enable easy comparison with the blue curve. 

In the work of Ding and Spelt |19j . phase field and level set models of a spreading droplet 
are compared for a range of Ohnesorge numbers (7.1 x 1CT 3 < Oh < 2.8 x 10 _1 ) making 
it a useful study for comparing with our work. The black dotted curve in figure [5] is a 
digitized curve of the lowest value of Oh simulated in Ding and Spelt. This particular 
simulation is selected for display here as it manifests the most pronounced oscillations. 
They simulate droplets with an initial contact angle of 7r/3 rad and an equilibrium contact 
angle of n/18 rad using an effective dimensionless slip length of A = 0.01 (A = 0.05 in our 
work). Despite these differences, the overall motion of the droplets agrees well qualitatively 
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for droplets with similar Ohnesorge numbers, although triple line motion was not seen to 
reverse direction in their work. In figure [6j the contact angle experiences a hysteresis loop 
in a similar fashion to the work of Ding and Spelt, which is reproduced in the black dotted 
curve. 

It has been conjectured [T91I2T] that the value of A Re* controls whether or not the spread- 
ing becomes oscillatory. In the simulations presented here, A Re* varies between 1 and 10 for 
the lowest value of Oh, but this is harder to determine for experimental systems. Hydrody- 
namic analysis of experimental data results in a slip length that can vary substantially for 
different materials (typically between 1 nm and 100 nm [7]). Using these bounds, a typical 
millimeter sized metal drop results in 0.01 < A Re* < 1 assuming a spreading rate of 1 m/s 
(in this work the spreading rate is ~ 50 m/s). It is interesting to note that for values of 
A Re* < 0.1, no oscillatory motion was seen in the work by Ding and Spelt |19| . In Schiaffino 
and Sorin [18] it is experimentally determined that the transition between under-damped 
oscillations to over-damped decay (no oscillations) occurs as Oh increases above 1 x 10~ 2 . 
This is seemingly confirmed in figure [4] where the curve that corresponds to Oh = 5.7 x 10~ 3 
has multiple oscillations, while the curve for Oh = 5.7 x 10 _1 has no oscillations. 

B. Comparison with experiments 

In figures [7] and [8] the triple-line radial position results from the present work for 
Oh = 5.7 x 10 -3 are compared with experimental results from Saiz and Tomsia [11] 
and Protsenko et al. [9|. These experiments are conducted at a high temperature (1100 °C) 
and exhibit fast spreading, which is either absent or undocumented in many other reactive 
wetting experiments J2H]- In Saiz and Tomsia, the experimental results are for Au and Cu 
droplets with an initial radius of 1 mm spreading on Ni and Mo substrates, respectively, 
while in Protsenko et al. the experiments are for Cu droplets of a similar size spreading on 
Si substrates. The reasonable quantitative agreement between the experimental and simu- 
lation results in figures [7] and [8] (within m 20% for the Cu-Mo combination) suggests that 
the spreading in the experimental systems is predominantly inertial in nature [TTj . 

The Cu on Mo spreading in figure [7] indicates oscillatory behavior at the end of the 
inertial regime, although there are only a handful of data points supporting this claim. 
Also, since the period of any oscillations is likely to be ~ 2U, a much greater duration of 



12 



experimental data is required for confirmation. The dissolutive case (black solid curve) in 
figure [8] clearly demonstrates oscillations of a similar period, amplitude and duration to the 
simulation results presented here as well as a contact angle hysteresis (not shown). It should 
be noted that oscillatory spreading also occurs in other systems such as water droplets on 
glass [18] . 



C. Non-equilibrium interface energy analysis 

The driving force for spreading on a planar substrate is often characterized by the spread- 
ing coefficient given by, 

S cqu (*) = 7s T - (iT + 7f7 cos 6 (f )) (16) 
where the 7 equ are equilibrium values of the interface energies and 9 is the observed contact 



angle. The utility of Eq. (16) is clearly limited to circumstances where the interface energies 
remain close to their equilibrium values during spreading. A number of authors (TU1 [2HH2I] 
have suggested that this limitation may be overcome by replacing the equilibrium interface 



energies with their instantaneous values in Eq. (16). This yields a new spreading coefficient 



S (t) = 7 W (t) - (%i (t) + % (t) cos 6 (t)) (17) 

where the 7 are instantaneous interface energies. In principal, the use of 7 rather than 7 equ 
provides a more accurate description of the driving force for spreading, particularly in the 
case where the timescale for spreading, ti, is much faster than the interface equilibration 
timescale, t^iff- Since the solid-fluid interface remains planar over the time scales of interest 
in the simulations, using a horizontal force balance alone and ignoring the vertical imbalance 



when deriving Eq. (17) can be viewed as a reasonable assumption. An alternative expression 
to Eq. (17) can be derived if the solid-fluid interface is non-planar using a more general 
Neumann's triangle horizontal and vertical force balance. In the following discussion, the 
expression used to calculate the 7 is described and then S is used to analyze the influence 
of £ on the spreading dynamics. 

It is a substantial advantage of our approach that we are able to develop an explicit 
expression for the instantaneous interface energies, allowing us to test the utility of S as a 
metric for spreading. In order to calculate S using the results of the present calculations, 
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we begin with two equivalent expressions for the equilibrium energy of a planar interface: 
7 = Jl~ hT|V Pl | 2 + e 2 T|Vp 2 | 2 + e,T|V0| 2 ] dl 

= 2 r_~ if-r- ( Pl - P f ) - & ( P2 - P f )] di 



where the oo superscript represents the value in the far field, and all the fields have equi- 



librium profiles. The equivalence of the expressions in Eq. (18) can be demonstrated by 



first writing down the Euler-Lagrange equation derived from the free energy functional in 



Eq. (A3) with additional Lagrange multiplier terms for the conservation of both species and 
then integrating once. We now assert that a plausible measure of the instantaneous interface 
energy is 

7 (*) = ^[eiT|Vpi| 2 + e 2 T|Vp 2 | 2 + e T|V0| 2 ] dl (19) 
where I is a line segment that both intersects and is normal to the interface being measured 



with j t dl > 5. All fields in Eq. ( 19 ) are measured at time t. In general, the quantity 7 is 



a useful heuristic when the gradients are confined to the interface region. The numerical 



integration of Eq. (19) is conducted at a distance of 25 from the triple-line location per- 
pendicular to each local interface over a distance of 1.55. The integration points on the 
respective interfaces are chosen to be as near to the triple-line location as possible while 
avoiding the large variations in the value of 7 that occur close to the triple-line location [I] . 
Clearly, we could have defined another instantaneous interface energy as, 

T *(t) = 2 J [f - r - {pi - PT) - 1% (P2 - p?)) dl (20) 

As one approaches equilibrium 7* — > 7, but dynamically the quantities are different. It 
would appear that 7* is less useful than 7, as 7* requires the fields to be near the far 
field (equilibrium) values at the integration limits extremes for the value to "make sense" 
as an interface excess quantity. It is instructive to observe the 7 behavior over time (see 
figure [9]). The values of 7 differ substantially from their equilibrium values for most of the 
simulation. The ^ sv appear independent of £, which is a reasonable expectation, as £ sets 
the liquid concentration. Increasing £ results in an increase in both 7^ and^. In figure [9] 
large oscillations can be observed in the solid-liquid interface energy (red curve). These 
oscillations are due to the spatially varying values of 7^ along the solid-liquid interface in 
conjunction with the oscillations in the j s i integration line location moving in unison with 
the triple-line location during the oscillatory phase of motion. 
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Using our definition of 7 and the apparent contact angle, 9, we can now calculate dynamic 



values of both S and >S equ , which are presented in figure 10 The curves decrease rapidly from 
their maximum value and become negative at about t = and then oscillate in conjunction 
with the triple-line radial position oscillations. Eventually, the values of S become quite 
small (< 10% of its original value for £ = 0) although the drop is still spreading. Assuming 
S quantifies the driving force for spreading, then the differences in S that occur for different 
values of £ at early times may explain both the deviations observed in the spreading extent 
during the inertial regime [t < tj) and the deviations in the oscillation amplitudes in figure [3j 
The small values of S when compared with S eqn at late times suggest that the spreading has 
become quasi-static in nature and is bound to the evolving values of the 7. The evolution 
of the 7 S / occurs on a time scale associated with tdm while the hydrodynamic adjustment 
of the contact angle occurs on a time scale associated with £j. Thus, the contact angle can 
adjust rapidly to balance the horizontal forces and suggests that the spreading is limited by 
interface equilibration at late times. 



D. Dissipation analysis 

Much of the literature surrounding droplet spreading is concerned with characterizing 
dissipation mechanisms from the point of view of an irreversible thermodynamic process [3, 
1321 133] . In this spirit, this section provides an analysis of the entropy production, yielding 
the magnitudes of the various dissipation mechanisms in our model, which should, in turn, 
provide guidance on the formulation of simplified models. The expression used here for the 
total entropy production rate is given by [34J , 

5prod = ^ W° - /£ c ) | 2 + ^ - + £ (d lUk + d kUl ) d lUk (21) 

where each term in the sum is a distinct dissipation mechanism (diffusion, solid interface 
relaxation, and viscous flow). 

The comprehensive overview of wetting by de Gennes [32] identified three main mecha- 
nisms for dissipation in spreading droplets: a viscous dissipation concerned with the "rolling 
motion" of the fluid within 100 /im of the triple line, a viscous dissipation in the precursor 
film and a highly localized dissipation at the triple line associated with "triple-line friction" . 
In the present work, the precursor film is absent, however, both viscous dissipation in the 
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bulk fluid and local triple-line dissipation are present, but are conflated within the viscous 



dissipation term in Eq. (21). In most models of droplets spreading, the chosen model for 
slip relaxation at the triple line influences the underlying dissipation mechanism for the 
spreading droplet. For example, a molecular kinetics model of slip generally implies a lo- 
cal triple-line dissipation, while a hydrodynamic model of slip, such as Cox's model [35] 
or Tanner's law [36], both examples of de Gennes' "rolling motion", implies non-localized 
dissipation jT], [33] . We are reminded that this model employs diffuse interfaces, and thus no 
explicit slip condition is postulated, but such slip is a direct consequence of the model. 



Figure 11 presents color contour plots of the entropy production rates at various times. 
The plots show the magnitude, location and mechanism of entropy production for the non- 
dissolutive case (the dissolutive cases are only slightly different). The color mapping is 

based on the max fSpRODj value for each image. For example, the 



rescaled in figure 
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total entropy production rate in figure 11 (d) is only 0.4% of the value in figure 11 (a). If we 



were considering a non-isothermal system, there would be a further term in expression 21 
containing temperature gradients, an effect not considered in this work. 

At very early times (not shown), the entropy production is highly localized at the solid- 
fluid interface region as </> locally equilibrates. Subsequently (not shown), pressure waves are 
observed as the liquid-vapor interface equilibrates, and viscosity is the dominant mode of 
dissipation. By t = O.ltj, the pressure waves have mostly subsided and the spreading is well 
under way. At this stage, the dominant dissipation mechanism remains viscous but is now 



highly localized at the triple-line. As the inertial time scale is approached in figure 11 (b), 
the dominant mechanism alternates between diffusive and viscous as the droplet oscillates 
during the U < t < 10ij stage. The viscous dissipation remains highly localized at the 
triple line, while the diffusive dissipation mostly occurs in the solid-liquid interface with 
some occurring along the solid- vapor interface. This correlates with figure [9j which shows 
that the solid-liquid interface is far from local equilibrium until much later times. At later 



times (figure 11 (c)), dissipation is mainly due to local interface equilibration along the solid- 



liquid and solid-vapor interface regions. The proportion of the numerically integrated value 



of J Sprodg^ for each term in Eq. (21) (diffusive, phase field, viscous) is (a) (0.51, 0.06, 
0.43), (b) (0.73, 0.04, 0.23), (c) (0.84, 0.02, 0.14) and (d) (0.85, 0.07, 0.08) for each subplot 



in figure 11 These proportions demonstrate the growing influence of diffusive dissipation 
and the reduction in viscous dissipation as the system transitions from the inertial regime 
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to the diffusive regime. 



E. Remarks 

The temporal adjustment to the equilibrium interface profiles is extremely complex and 
intimately related to the interface width and the interpolated values of the dynamic coef- 
ficients {y and M), see Eqs. (|5]) and (|6]). The choice for the interpolation parameter a in 
Eqs. ([5]) and ^ biases the coefficients to have values close to the bulk fluid values in the 
interface region facilitating the fastest interface dynamics possible within the bounds set by 
the bulk values. The parameter a is tuned to a value of 4, as larger values do not increase 
the interface equilibration rate while smaller values considerably reduce the equilibration 
rate. 

The equilibration of the density and phase field interface profiles is fast compared to that 
of the concentration field. The interface profile of the density field, p, is adjusted rapidly 
by hydrodynamics alone, while the interface profile of the concentration field, p\j p, requires 
inter-diffusion between the bulk phases and the interface regions. This compositional relax- 
ation could, in principle, be as slow as the diffusion time scale, tdiff (see Table |ll|) , although 
the connection is imprecise, as this quantity is associated primarily with the motion of the 
interface due to dissolution (melting) rather than the relaxation of compositional profiles 
within the interface. The solid interfaces equilibrate slowly, compared to the liquid-vapor 
interface, as seen in figure [9j We expect that the observed interface relaxation time is unreal- 
istic, when compared with experimental studies of metallic systems, as our chosen interface 
width of 5 = 100 nm is much larger than the 8 ~ 1 nm typical of metals. This is a shortcom- 
ing of this treatment, and results in an unphysical time scale for local interface equilibration. 
Further analysis of the relationship between S and the equilibration rate is required, though 
this analysis is beyond the scope of this work. The limitation of requiring 5/R ~ 0.1 im- 
posed by the available computation resources does not detract from the analysis presented 
in this section with respect to the reduced spreading when £ is increased, the qualitative 
description of the spreading regimes and oscillations, or the quantitative comparisons with 
experiments. 
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VI. CONCLUSION 



This paper presents results from a model of dissolutive spreading simulated in a parame- 
ter regime where inertial effects are initially dominant. The triple-line motion demonstrates 
good agreement with the 0(t~^ 2 ) inertial spreading rate at early times. The model also 
generates oscillations characteristic of the transition from inertial to viscous or diffusive 
spreading. Subsequent analysis indicates that a force balance involving the instantaneous 



interface energies evaluated using the expression in Eq. (19) can explain the variation in 
spreading between the hydrodynamic and dissolutive cases. At late times, after inertial 
effects have ceased, the contact angle derived from the instantaneous interface energies is 
within 0.005 rad of the measured contact angle suggesting that the local interface equilibra- 
tion mechanism is controlling the spreading. Analysis of the dissipation mechanism via the 
entropy production expression demonstrates that dissipation occurs at the triple line during 
the inertial stage, but transitions to the solid-fluid interfaces during the oscillatory stage 
consistent with the instantaneous interface energy analysis. Overall, the simulation results 
show good quantitative and qualitative agreement with a number of experimental results 
when time is scaled with the inertial time scale. 

Modeling droplets that have both a realistic interface width and include inertial effects is 
impractical with current computational resources (at least for the model presented herein) 
and may require years of real time computation on large parallel clusters. In this work, to 
reduce the required compute time, the use of a realistic interface width has been sacrificed in 
order to preserve the inertial effects. This has the consequence of increasing the simulation 
time required for the local equilibration process across the solid-fluid interface as discussed 



in section |V C| Although, this process has a longer duration than physically appropriate in 
the present work, a time regime over which the controlling mechanism for spreading is the 
local interface equilibration may be entirely physical. It is noted in Protsenko et al. [U] that 
the diffusive stage may occur in two separate parts. The first part is surmised to be the 
solid-liquid interface equilibration process and takes approximately an order of magnitude 
longer than the inertial time scale, which is faster than occurs here, but very similar in 
nature. The second part is the melting of the substrate, which is included in this model, 
but not observed as it occurs over a time scale longer than the total duration of a typical 
simulation. 
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Further work may involve both direct comparison with molecular kinetics theory and 
more detailed analysis of the impact of the interface width on the spreading dynamics. 
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Appendix A: Derivation of the Governing Equations 



In this section the underlying thermodynamic and constitutive relationships required for 
the derivation of Eqs. §2§, ^ and Q are presented. 

As previously outlined, the fluid phases are represented by a binary, van der Waals 
equation of state and the solid phase is represented by a simple linear compressive and 
tensile equation of state that ignores all shear stress. The van der Waals equation of state 
is given by, 

' 2 

n 



P - (e 1 X 1 + e 2 X 2 )j (V - vn) = nRT (Al) 
where X± and X 2 are the concentrations of each component, n is the number of moles and 



V/n = m/p. All other parameters used in Eq. (Al) are defined in section |Hj Eq. (Al) can 
be related to the ideal gas law, but has modified pressure and volume terms to account for 
the long range attraction of molecules and volume exclusion, respectively j25j|37j. The solid 
equation of state is given by, 



PV„ = 2Bn 



Vs-V 
V. 



(A2) 



where V s /n = m/p r s ei . The Helmholtz free energies given in Eqs. (|7j) and ^ are derived 



from (Al ) and (\A2y, respectively, using the thermodynamic identities given in Eqs. (10), (11 ) 



and ([9]). In order to derive Eqs. (]2j, ^ and Q, it is necessary to postulate a form for the 
free energy functional. 

As in reference |38j, standard non-classical diffuse interface expressions for p±, p 2 and 
are used, which results in a functional of the form, 



dV 



2 1 2 

Using standard dissipation arguments [38], Eqs. ^ and ^ are derived using, 

5F 



(A3) 



at 



-Mm 



8(f) 
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and 

dp 
~dt 



- + dj (iijpx) = -djhj 



and similarly for component 2. The fluxes are given by, 



Jij = -J23 = —Mdj 



T 



where, 

l J i 



Spi 
and 

nc _ *F 

The form of the stress tensor required to derive the momentum equation is given by, 

<Jij = v (djm + diUj) + Uj 

using the standard assumption that the bulk viscosity, A, is related to the shear viscosity via 
A = The tensor, tij, is derived from a conservation law (djtij = 0) based on Noether's 

theorem [39J. The expression for is given by, 

f « = 9 ^-^afe) (A4 > 

where 

g NC = f NC + Pi\i + P2\2 (A5) 

The non-classical Gibbs free energy, g NC , is the form of the free energy that includes La- 
grange multipliers for conservation of species 1 and 2. The Lagrange multipliers for each 
species are equal to Ai = —fi^ c and A2 = —^2° * n equilibrium using the variational deriva- 



tive of J g NC dV with respect to p\ and pi- Using Eqs. (A4) and (A5) the form for d{tij used 
in Eq. Q can be derived, 

djUj = -Pid^ c - ptditf ° - (A6) 
Appendix B: Numerical Approach 

In general, even for compressible systems, many conventional algorithms use the pressure 
field as the independent variable rather than the density field. This approach is thought 
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to have more robust convergence properties |40j at low Mach numbers due to the weak 
dependence of pressure gradients on density, but the convergence properties deteriorate at 
higher Mach numbers. In this work, due to the non-trivial nature of the pressure-density 
relationship, an inversion of this relationship would be impractical and it is more natural to 
solve for the density field rather than the pressure field. Due to the mesh collocation of the 
density and velocity fields, an interpolation scheme, known as Rhie-Chow interpolation jjl] , 
is employed to ensure adequate velocity-pressure coupling. 

The calculation of triple-line velocities is necessarily noisy, with fluctuations on a timescale 
of Ax/U, where Aa; is the fixed grid spacing. In figure [5j the curves are constructed using a 
20 point boxcar (equally-weighted) averaging scheme collected at every 10 time steps during 
the simulation. We note that the sign changes in the blue curve (U < t < 10U) in figure [5] 
correspond to the triple-line oscillations, and are not due to the averaging scheme. The 
velocity fluctuations will be small when U is large. Indeed, at early times, when t/ti < 1, U 
is relatively large and the results are smooth. At later times, when t/ti > 10, the averaging 
scheme does not smooth out the noise, as the spreading rate is greatly reduced. This can 
be seen in the noisy behavior at long times for the Oh = 5.7 x 10 -3 curve (blue) in figure [5} 
The noise in the low velocity regime of figure [6] also reflects this behavior. 

The measurements for 9 are calculated using the tangent to the liquid-vapor interface 
at a distance of 1.3 5 from the triple-line location. In general, this distance results in a 
reasonable approximation to the apparent contact angle. 

1. Parasitic Currents 

Parasitic currents are a common source of numerical errors when computing flows with 
interface energy driving forces that have large Ca. Typically, for the systems of interest in 
this paper, Ca ~ 10 -2 , but parasitic velocities were still found to be a source of numerical 
error, particularly when trying to evaluate equilibrium solutions. Parasitic currents are 
characterized by quasi-steady flow fields that do not dissipate over time despite the system 
reaching equilibrium in all other respects. This can result in equilibrium errors in both 
the density and concentration fields. Jamet et al. |42J as well as other researchers have 
demonstrated that parasitic currents can be eliminated by recasting the momentum equation 
in a form that only conserves momentum to the truncation error of the discretization rather 
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than machine precision. The form of the momentum equation that eliminates parasitic 
currents is written in terms of the chemical potentials and is given by, 

d(pUi 



dt 



+ dj (pUiUj) = dj [djUi + diUj}) - pid^i ° - p 2 dip 2 



,NC 



(Bl) 



for binary liquid- vapor system. The discretized form of Eq. (Bl ) is known as an energy con- 
serving discretization in contrast to the momentum conserving discretization, which results 



when the momentum equation is written in terms of the pressure (see Eq. (13)). 



2. Convergence 

Some simulations in this paper are tested for convergence with grid sizes of 180x125, 
360x250 and 720x500 using the triple-line and drop height positions against time as the met- 
rics for convergence. Production runs for the results presented use 360x250 grids. Details 
of these convergence tests can be found in [23J . Convergence at the n th time-step is achieved 
when the k th iteration within the time step satisfies the residual condition < 1 x 10 _1 

for each of the equations where is the L 2 -norm of the residual at the the k th iteration of 
the n th time step. Further decreases in the residual make little difference to the dynamic po- 
sitions of the drop height and triple-line. Numerical calculations indicate that, in the course 
of a simulation, Ma ranges from values that require compressible flow solvers (density based 
with Ma > 2 x 10 _1 ) to values for which compressible flow solvers have trouble with accuracy 
and convergence for traditional segregated solvers (Ma < 2 x 10 _1 ). The shift to low Ma 
generally occurs when the system is quite close to equilibrium and is not believed to affect 
the dynamic aspects of the simulation, which are of most interest in this paper. In general, 
for low Mach number flows, preconditioners are used to improve the convergence properties 
of segregated solvers. In this work, it was found that using a coupled solver along with a 
suitable preconditioner greatly improved the convergence properties. The preconditioners 
are available as part of the Trilinos software suite [13]. The coupled convergence properties 
can be further improved by employing physics based preconditioners that change the nature 
of the equations based on the value of Ma [44J , but are not used in this work. 
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Parameter 


Value 


Unit 


V S 


2.0xl(T 3 


kg/ (s • m) 




2.0xl0 4 


kg/ (s • m) 


ei 


2.0xl0- 16 


m 7 /(K-kg-s 2 ) 




2.0xl0~ 16 


m 7 /(K-kg-s 2 ) 


T 


6.5xl0 2 


K 


m 


1.18X10" 1 


kg/mol 


R 


8.31 


J/(K-mol) 


V a 


1.0 




ei 


-4.56X10" 1 


J • m 3 /mol 2 




-4.56X10- 1 


J • m 3 /mol 2 


V 


1.3xl0~ 5 


m 3 /mol 


A 1 


2.83xl0 4 


J/mol 


A 2 


5.64xl0 4 


J/mol 


pf 


7.84xl0~ 5 


kg/m 3 


B 


2.02xl0 5 


J/mol 


W 


1.27xl0 5 


N/m 2 


e 4> 


l.OxlO" 9 


N/K 




l.OxlO 4 


K-m 2 /(N-s) 


M f 


1.0xl0~ 7 


kg • s • K/m 3 


M s 


l.OxlO- 11 


kg • s • K/m 3 


Rq 


l.OxlO" 6 


m 


5 


l.OxlO" 7 


m 


equ 
Pi 


7.35 xlO 3 


kg/m 3 



TABLE I: Various parameter values. 
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Time scale 


Symbol Expression 


Value (s) 


capillary 


tc 


VfRo/jlv 


1.05xl(T 10 


phase field 




S 2 /e (j) M (f> 


l.OxlO^ 9 


inertial 


U 


a / P? qU -Ro/7/jj 


1.97xl0- 8 


convection 


ta 


R /U 


1.97xl(T 8 


viscous 


tu 




3.68xl0~ 6 


interface diffusion 


tdiS 


5 2 /4K 2 D f 


7.69xl0" 4 


bulk diffusion 


t d 


Rl/Dj 


1.04xl0~ 3 


solid deformation 


ts 




1.05xl0~ 2 


instantaneous convection t* 


Ro/U* 





TABLE II: Complete list of time scales referred to in this paper. 



Parameter 


Symbol Expression 


Value 


Peclet number 


Pe 


URo/Df = t d /t a 


5.31 x 10 4 


Reynolds number 


Re 


URop^/i/f = t v /t a 


1.87 x 10 2 


Weber number 


We 


Re Ca = t v t c jt\ 


1.0 


Mach number 


Ma 


U/c 


5.72 x 1(T 2 


Unnamed 


Q 


m^/RTp^Ro 


5.64 x 1(T 2 


effective dimensionless slip length A 


S/Ro 


5.0 x 10~ 2 


Capillary number 


Ca 


U Vfhlv = tc/ta 


5.35 x 10~ 3 


Ohnesorge number 


Oh 


V / C&/Re = t c /t i 


5.35 x 10-3 


instantaneous Reynolds number 


Re* 


U*R 0P c r/vf = U/t* a 




instantaneous Capillary number 


Ca* 


U*VfHv = tc/t* a 





TABLE III: Relevant dimensionless numbers. 
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FIG. 1: The phase diagram for the system of parameters presented in Table [IJ Each region rep- 
resents a possible equilibrium state for a mixture of solid (S), liquid (L) and vapor (V) phases. 
The red dots represent the initial conditions for the £ = 0.1 simulation discussed in section [TV] 
The black dot marks the liquid equilibrium condition. The liquid and vapor phases are thick in 
component 1 while the solid phase is thick in component 2. 
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r (Aim) 



FIG. 2: Sequential configurations of the liquid- vapor and solid-fluid interfaces for £ = and 
Oh = 5.7 x 10 -3 with darker tones indicating later times. The curves demonstrate the extreme 
inertial effects on the droplet. The droplet starts as a sphere in tangent contact with the substrate. 
The drop height then rises considerably as the capillary wave initiated from the triple line arrives 
at the top of the droplet. Although large amplitude (« Rq/5) oscillations occur in the triple- line 
position, the largest contact angle oscillation is only « 0.037T rad. 
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FIG. 3: The spreading radius versus time for various values of £ with Oh = 5.7 x 10 3 . As £ 
increases, the spreading rate and extent of spreading is slightly reduced. 
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FIG. 4: The spreading radius, r#, versus time for various values of Oh with £ = 0. The oscillations 
are eliminated for the largest value of Oh. 
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FIG. 5: The dimensionless spreading rate against the dimensionless time with varying Oh and 
£ = 0. The spreading occurs in three distinct intervals.. The sign changes in the blue curve 
correspond to the triple-line oscillations during the transition from the inertial to the diffusive 
regime. 
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FIG. 6: The observed contact angle against Ca* for £ = and Oh = 5.7 x 10 
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FIG. 7: The radial position of the triple line scaled against the final radial position, Rf, against 
time (scaled with U) for £ = 0, Au-Ni experimental results and Cu-Mo experimental results. The 
experimental results are digitized from Saiz et al. [7J [11]. The inertial time scale, t% for the Au-Ni 
system is calculated using p = 1.1 x 10 4 kg/m 3 , 7 = 1.0 J/m 2 , R = 1 x 1CT 3 m. The inertial 
time scale for the Cu-Mo system is calculated using p = 8.9 x 10 3 kg/m 3 , 7 = 1.3 J/m 2 and 
Rq = 1 x 1CT 3 m. The value of ij is 1.9 x 1CT 8 s for this work, 3.4 x 1CT 3 s for the Au-Ni system 
and 2.6 x 10 -3 s for the Cu-Mo system. This figure shows the reasonable agreement between the 
simulation and experimental data when scaled by the inertial time scale and the agreement with 
the {t/Ati) 1 / 2 spreading rate. 
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FIG. 8: The spreading radius versus time for £ = and £ = 0.9. The black curves are Cu-Si 
experiments digitized from Protosenko et al. The inertial time scale, U, for the Cu-Si system is 
calculated using p = 8.9 x 10 3 kg/m 3 , 7 = 1.3 J/m 2 and Rq = 8.2 x 1CT 4 m. The value of t% is 
1.9 x 10 -8 s for this work and 2.0 x 10~ 3 s for the Cu-Si system 
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FIG. 9: The instantaneous interface energies 7 plotted against time for £ = and £ = 0.9. Both 
7^ and j s i are larger for the £ = 0.9 curve. 
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FIG. 11: Contour plots of the entropy production rate at (a) t = O.ltj, (b) t = ti, (c) t = 10tj 



and (d) t = 20U. The color intensity represents the magnitude of either y/ Sprod (less focused) 
on the left panel or \J Sprod (more focused) on the right panel. The colors represent the specific 



entropy production mechanism given by the terms in Eq., (21 ) (diffusive, phase field, viscous), with 
red, green and blue representing the first (diffusion), second (solid interface relaxation) and third 
(viscous flow) terms, respectively. 
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